subroutine initial()
use type
use var
implicit none
integer :: i
logical :: flag = .true.

  bnodes=0 
  boundary = .false.
  drop = .false.
  lb2index = 0
  call grid()
  q = 0.0
  neighbors = 0
  call find_neighbors()
  call find_normal0()
  call anchoring()
  call find_boundary_neighbors0()
  call q_initial()
!  call qmga(0,0)

contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Subroutine grid
use type
use var
implicit none
integer(i4b) :: i,j,k,label
integer(i4b) :: cont_b = 0 , cont_d = 0, R
real(dp) :: checkx, checky, checkz, check1, check2, check3
real(dp) :: Lx, Ly, Lz
!********************************************
!loop for surface and inside nodes
center = int(dime/2) + 1
R = center(1)-2
cont_b = 0
cont_d = 0
  do k=1,dime(3)
     do j=1,dime(2)
        do i=1,dime(1)

           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
           if (channel) then 
              if (k==1.or.k==dime(3)) then
                  boundary(i,j,k) = .true.
                  cont_b = cont_b + 1
                  lb2index(label) = cont_b
                  if (k==1) boundary0(i,j,k) = .true.
              endif
           elseif (droplet) then
              Lx = i-center(1)
              Ly = j-center(2)
              Lz = k-center(3)

              checkx = (Lx)/real((R+0.5))
              checkx = checkx*checkx
              checky = (Ly)/real((R+0.5))
              checky = checky*checky
              checkz = (Lz)/real((R+0.5))
              checkz = checkz*checkz
              check1 = checkx+checky+checkz
           

              if (check1<= 1.d0) then
                 drop(i,j,k) = .true.
                 cont_d = cont_d + 1
              endif
            elseif (cylinder) then
              Lx = i-center(1)
              Ly = j-center(2)

              checkx = (Lx)/real((R+0.5))
              checkx = checkx*checkx
              checky = (Ly)/real((R+0.5))
              checky = checky*checky
              check1 = checkx+checky           

              if (check1<= 1.d0) then
                 drop(i,j,k) = .true.
                 cont_d = cont_d + 1
              endif


           endif  
        enddo
     enddo
  enddo

if (droplet) then
  do k=2,dime(3)-1
     do j=2,dime(2)-1
        do i=2,dime(1)-1

           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)

           if (drop(i,j,k).and.(drop(i+1,j,k)==.false..or.drop(i-1,j,k)==.false. &
                            .or.drop(i,j+1,k)==.false..or.drop(i,j-1,k)==.false. &
                            .or.drop(i,j,k+1)==.false..or.drop(i,j,k-1)==.false.)) then
              boundary(i,j,k) = .true.
              cont_d = cont_d - 1
              cont_b = cont_b + 1
              lb2index(label) = cont_b
             
           endif

        enddo
    enddo
 enddo
endif

if (cylinder) then
  do k=1,dime(3)
     do j=2,dime(2)-1
        do i=2,dime(1)-1

           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)

           if (drop(i,j,k).and.(drop(i+1,j,k)==.false..or.drop(i-1,j,k)==.false. &
                            .or.drop(i,j+1,k)==.false..or.drop(i,j-1,k)==.false.)) then
              boundary(i,j,k) = .true.
              cont_d = cont_d - 1
              cont_b = cont_b + 1
              lb2index(label) = cont_b
             
           endif

        enddo
    enddo
 enddo
endif

if (droplet.or.cylinder) then
 do k=1,dime(3)
     do j=1,dime(2)
        do i=1,dime(1)

           if (boundary(i,j,k)) drop(i,j,k) = .false.

        enddo
    enddo
 enddo
endif

bnodes = cont_b
dnodes = cont_d

end subroutine grid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Subroutine find_neighbors
use type
use var
implicit none
integer(i4b) :: i, j, k, ii 
integer :: x, y, left, right, label_l, label_r, label
integer, dimension(3) :: pos, pos_l, pos_r

  do k=1,dime(3)
     do j=1,dime(2)
        do i=1,dime(1)
           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
           pos(1) = i ; pos(2) = j ; pos(3) = k
           Do x = 1 , 3
              left = pos(x) - 1
              right = pos(x) + 1
              if (left<1) left = dime(x)
              if (right>dime(x)) right = 1
              Do y = 1 , 3
                 if (y/=x) then
                    pos_l(y) = pos(y)
                    pos_r(y) = pos(y)
                 else
                    pos_l(y) = left
                    pos_r(y) = right
                 endif
              enddo
              label_l = 1 + (pos_l(1)-1) + (pos_l(2)-1)*dime(1) + (pos_l(3)-1)*dime(1)*dime(2)
              label_r = 1 + (pos_r(1)-1) + (pos_r(2)-1)*dime(1) + (pos_r(3)-1)*dime(1)*dime(2)
              ii = (x-1)*2 + 1
              neighbors(ii,label) = label_l
              neighbors(ii+1,label) = label_r
           enddo
        enddo
     enddo
  enddo
end subroutine find_neighbors
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine find_normal0()
use var
use type
implicit none
integer :: i,j,k, indx, label
real(dp) :: x,y,z
real(dp) :: theta, phi, r
if (allocated(normal0)) deallocate(normal0)
if (allocated(tangent0)) deallocate(tangent0)
allocate(normal0(3,bnodes),tangent0(3,bnodes)) 
normal0=0.d0 ; tangent0=0.d0; indx = 1

if (channel)  then
   do k=1,dime(3)
      do j=1,dime(2)
         do i=1,dime(1)
            if (boundary(i,j,k)) then
               label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
               indx = lb2index(label)
               if (k==1) normal0(3,indx) = 1.d0
               if (k==dime(3)) normal0(3,indx) = -1.d0
            endif
         enddo
      enddo
   enddo
endif

if (droplet) then
   do k=1,dime(3)
      do j=1,dime(2)
         do i=1,dime(1)
            if (boundary(i,j,k)) then
               x = real(i - center(1)); y = real(j - center(2)); z = real(k - center(3))
               r = sqrt(x*x + y*y + z*z)
               theta = acos(z/r)
               if (x/=0) then
                  phi = atan(y/x)
               else
                  phi = Pi /2.d0
               endif
               if (x<0.d0.and.y<0.d0) then
                  phi = phi + Pi
               elseif (x<0.d0.and.y>0.d0) then 
                  phi = phi + Pi
               elseif (x>0.d0.and.y<0.d0) then
                  phi = phi + 2*Pi
               endif
               if (phi==0.d0.and.x<0) then
                  theta = Pi - theta
               endif
               if (phi==Pi/2.d0.and.y<0) then
                  theta = Pi - theta
               endif
               normal0(1,indx) = sin(theta)*cos(phi)
               normal0(2,indx) = sin(theta)*sin(phi)
               normal0(3,indx) = cos(theta)
               tangent0(1,indx) = cos(theta)*cos(phi)
               tangent0(2,indx) = cos(theta)*sin(phi)
               tangent0(3,indx) = -sin(theta)
               if (phi==0.d0.and.x<0) then
                  normal0(:,indx) = -normal0(:,indx)
               endif
               if (phi==Pi/2.d0.and.y<0) then
                  normal0(:,indx) = -normal0(:,indx)
               endif
               indx = indx + 1
            endif
         enddo
      enddo
   enddo
   normal0 = - normal0 ! the unit vector of droplet inside 
endif
if (cylinder) then
   do k=1,dime(3)
      do j=1,dime(2)
         do i=1,dime(1)
            if (boundary(i,j,k)) then
               x = real(i - center(1)); y = real(j - center(2))
               r = sqrt(x*x + y*y)
               if (x/=0) then
                  phi = atan(y/x)
               else
                  phi = Pi /2.d0
               endif
               if (x<0.d0.and.y<0.d0) then
                  phi = phi + Pi
               elseif (x<0.d0.and.y>0.d0) then 
                  phi = phi + Pi
               elseif (x>0.d0.and.y<0.d0) then
                  phi = phi + 2*Pi
               endif
               normal0(1,indx) = cos(phi)
               normal0(2,indx) = sin(phi)
               normal0(3,indx) = 0.d0 
               tangent0(1,indx) = -sin(phi)
               tangent0(2,indx) = cos(phi)
               tangent0(3,indx) = 0.d0
               if (phi==0.d0.and.x<0) then
                  normal0(:,indx) = -normal0(:,indx)
               endif
               if (phi==Pi/2.d0.and.y<0) then
                  normal0(:,indx) = -normal0(:,indx)
               endif
               indx = indx + 1
            endif
         enddo
      enddo
   enddo
   normal0 = - normal0 ! the unit vector of droplet inside 
endif
end subroutine find_normal0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
subroutine anchoring
  use type
  use var
  implicit none
  integer(i4b) i,j,k, label, indx 
  real(dp) max_e, s, dr
  real(dp), dimension(3):: r, planar = 0.d0, eigenvalues, rand 
  
  if (allocated(qb0)) deallocate(qb0)
  allocate(qb0(6,bnodes))  ! boundary q for surface
  indx = 1

  do k=1,dime(3)
     do j=1,dime(2)
        do i=1,dime(1)

           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
           if (channel) then 
           if (boundary(i,j,k) == .true.) then
              
              !type of channel
              planar = 0
              if (k==1) then
                 if (Stype/="zz") planar(1) = 1
                 if (Stype=="zz") planar(3) = 1
              else
                 if (Stype=="xx") planar(1) = 1
                 if (Stype=="xy") planar(2) = 1
                 if (Stype=="xz") planar(3) = 1
                 if (Stype=="zz") planar(3) = 1
              endif

              call director_tensor( Sbulk, planar, qb0(1:5,indx))
              ! Checking all the values of the tensor  
              call eigen( qb0(1:5,indx), eigenvalues, r )   ! max_e is the maximun eigenvalue and r is the eigenvector
              s = maxval( eigenvalues )*3.d0/2.d0
              if ( abs(s - Sbulk) > 1.d-10 ) then
                 print *, '--> Problem with anchoring 1', i, r
                 stop
              end if
              dr = sqrt( dot_product( r,r ) )
              if ( abs(dr - 1.d0) > 1.d-10 ) then 
                 print *, '--> Problem with anchoring 2', i, r
                 stop
              end if
              
              qb0(6,indx) = - qb0(1,indx) - qb0(4,indx)
              q(:,label) = qb0(:,indx)
              indx = indx+1
           endif
           elseif (droplet.or.cylinder) then
           if (boundary(i,j,k) == .true.) then
              
              !type of channel
              planar = 0
              if (Stype=="xx") planar = tangent0(:,indx) !plannar
              if (Stype=="zz") planar = normal0(:,indx)  !homotropic
              if (Stype=="xy".and.cylinder) planar(3) = 1.0  !homotropic
               
              if (itype==3) then
              call RANDOM_NUMBER( rand )
              rand = rand - 0.5d0
              dr = sqrt( dot_product( rand, rand ) )
              rand = rand/dr
              planar = rand
              endif

              call director_tensor( Sbulk, planar, qb0(1:5,indx))
              ! Checking all the values of the tensor  
              call eigen( qb0(1:5,indx), eigenvalues, r )   ! max_e is the maximun eigenvalue and r is the eigenvector
              s = maxval( eigenvalues )*3.d0/2.d0
              if ( abs(s - Sbulk) > 1.d-10 ) then
                 print *, '--> Problem with anchoring 1', i, r
                 stop
              end if
              dr = sqrt( dot_product( r,r ) )
              if ( abs(dr - 1.d0) > 1.d-10 ) then 
                 print *, '--> Problem with anchoring 2', i, r
                 stop
              end if
              
              qb0(6,indx) = - qb0(1,indx) - qb0(4,indx)
              q(:,label) = qb0(:,indx)
              indx = indx+1
           endif
        endif
       
        enddo
     enddo
  enddo
  

end subroutine anchoring
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine q_initial
  use type
  use var
  implicit none
  integer(i4b) i,j,k, label
  real(dp) max_e, s, dr
  real(dp), dimension(3):: r = 0.d0, rand = 0.d0, eigenvalues
  real(dp), dimension(5) :: q_temp 
  do k=1,dime(3)
     do j=1,dime(2)
        do i=1,dime(1)

           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
           
           if (boundary(i,j,k) == .false. ) then  
              call RANDOM_NUMBER( rand )
              rand = rand - 0.5d0
              dr = sqrt( dot_product( rand, rand ) )
              rand = rand/dr
              if (itype == 2) then
                 rand = 0.0
                 if (channel) rand(1) = 1.0
                 if (droplet) rand(3) = 1.0
              endif
              call director_tensor( Sinit, rand, q_temp )
              q(1:5,label) = q_temp
              ! Checking all the values of the tensor  
              call eigen( q_temp, eigenvalues, r )   ! max_e is the maximun eigenvalue and r is the eigenvector
              s = maxval( eigenvalues )*3.d0/2.d0
              if ( abs(s - Sinit) > 1.d-10 ) then
                 print *, '--> Problem with anchoring 1', i, r
                 stop
              end if
              dr = sqrt( dot_product( r,r ) )
              if ( abs(dr - 1.d0) > 1.d-10 ) then 
                 print *, '--> Problem with anchoring 2', i, r
                 stop
              end if

           endif
        
        enddo
     enddo
  enddo
end subroutine q_initial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine find_boundary_neighbors0()
use var
use type
implicit none
integer :: i,j,k, indx, label
integer :: deltax,deltay,deltaz

if (allocated(boundary_neighbors0)) deallocate(boundary_neighbors0)
allocate(boundary_neighbors0(6,bnodes))
boundary_neighbors0 = 0
indx = 1

do k=1,dime(3)
   do j=1,dime(2)
      do i=1,dime(1)
         label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2) 
         if (channel) then
         if (k==1) then
            indx = lb2index(label)
            boundary_neighbors0(5,indx) = 1 + (i-1) + (j-1)*dime(1) + (k-1+1)*dime(1)*dime(2)
            boundary_neighbors0(6,indx) = 1 + (i-1) + (j-1)*dime(1) + (k-1+2)*dime(1)*dime(2)
         elseif (k==dime(3)) then
            indx = lb2index(label)
            boundary_neighbors0(5,indx) = 1 + (i-1) + (j-1)*dime(1) + (k-1-1)*dime(1)*dime(2)
            boundary_neighbors0(6,indx) = 1 + (i-1) + (j-1)*dime(1) + (k-1-2)*dime(1)*dime(2)
         endif
         elseif (droplet.or.cylinder) then
            if (boundary(i,j,k)==.true.) then
            indx = lb2index(label)
            if (normal0(1,indx)>=0) deltax=1
            if (normal0(1,indx)< 0) deltax=-1
            if (normal0(2,indx)>=0) deltay=1
            if (normal0(2,indx)< 0) deltay=-1
            if (normal0(3,indx)>=0) deltaz=1
            if (normal0(3,indx)< 0) deltaz=-1
            boundary_neighbors0(1,indx) = 1 + (i-1+deltax) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
            boundary_neighbors0(2,indx) = 1 + (i-1+2*deltax) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
            boundary_neighbors0(3,indx) = 1 + (i-1) + (j-1+deltay)*dime(1) + (k-1)*dime(1)*dime(2)
            boundary_neighbors0(4,indx) = 1 + (i-1) + (j-1+2*deltay)*dime(1) + (k-1)*dime(1)*dime(2)
            boundary_neighbors0(5,indx) = 1 + (i-1) + (j-1)*dime(1) + (k-1+deltaz)*dime(1)*dime(2)
            boundary_neighbors0(6,indx) = 1 + (i-1) + (j-1)*dime(1) + (k-1+2*deltaz)*dime(1)*dime(2)
         endif
         endif
     enddo
   enddo
enddo
end subroutine find_boundary_neighbors0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
end subroutine initial
